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ABSTRACT 


This  report  documents  the  computer  code  WAVEAMP,  a  program  to  compute  the 
Kelvin  wave  amplitudes  using  the  Nuenann-Kelvin  approxi-ation.  To  run  WAVEAMP 
the  source  strength  on  the  hull  surface  must  first  be  determined  using  a 
Neumann-Kelvin  code. 

Example  output  is  presented  for  the  Quapaw  hull  at  three  Froude  numbers. 
The  calculations  include  both  near-field  and  far-field  waves,  and  longitudinal 
wave  cuts. 
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INTRODUCTION  AND  OUTLINE 


This  report  is  divided  into  four  sections.  A  brief  theoretical 
background  and  problem  formulation  is  included  in  Section  I.  The  numerical 
method  used  is  described  in  Section  II.  The  method  results  in  a  considerable 
reduction  in  computer  time  at  the  expense  of  slight  loss  of  accuracy.  In 
Section  III  results  of  predicted  wave  amplitude  for  the  QUAPAW  hull  are 
presented.  Finally,  user's  instructions  for  program  WAVEAMP  are  described  in 
Section  IV.  A  program  listing  is  included  in  the  Appendix. 
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I.  THEORETICAL  BACKGROUND 


In  this  section,  a  brief  outline  of  the  theoretical  background  needed  for 
the  computation  of  ship  waves  is  presented.  Details  of  the  material  included 
here  can  be  found  in  (1]  and  [21  • 

1 .1  Problem  Formulation 

In  the  following,  a  right  handed  coordinate  system  is  used  with  the  origin 
at  the  intersection  of  the  calm  water  plane  and  midship.  The  positive  x-axis 
points  out  the  bow  and  the  z-axis  in  position  upward. 

The  total  potential  for  a  thick  ship  advancing  with  a  constant  speed  V 
in  an  otherwise  cairn,  inviscid  and  incompressible  fluid  is  given  by 

♦  ■  -  VX  +  <1>  t  (  1  ) 

where  <*  is  the  perturbation  potential  which  satisfies  Laplace's  equation 

=0  (2) 
The  linearized  free-surface  boundary  condition  is 

*xx  +  k0>2  =  0  on  z  =  0  ,  (3) 

where  kQ  =  g/V°  .  A  fundamental  solution  of  (2)  and  (3)  is  the  Green 
function,  given  by 

1  1  kQ  ~  «•>  exp^k  [z+'+i  (x-i  )cost<+i  (y-n  Isin^1 

G(P,Q)  - - +  —  Lin  r  d°  r  dk - - -  (4) 

r  r'  i  i'-»0  -ii  0  k0  -  k  cos’0  -  in  cosc 

where  P  =  (x,y,z)  is  the  field  point,  Q  =  is  a  point  source  of 

strength  -4m  ,  and 

r,  r’  =  [  { x-r )-  +  (y-n )?  +  (z  =  .  (5) 

In  the  usual  method  of  potential  theory.  Green's  theorem  can  be  applied  to 
a  large  volume  of  fluid  containing  the  body,  and  extending  to  infinity  both 
laterally  and  in  depth.  The  following  result  for  the  perturbation  potential  in 
terms  of  a  pure  source  density  distribution  r  is  obtained: 


2- 
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1  1 

4>(P)  - - //  G(P,Q)o(Q)dS(Q) - /  vx  2( p, Q)g ( 2)dn  ,  (6 ) 

4*  SH  4”k0  CH 

where  vx  is  the  x-component  of  the  normal  vector  out  of  the  fluid  into  the 
source  domain,  SH  is  the  wetted  surface  of  the  body,  and  CH  is  the 
intersection  of  the  hull  and  the  calm  water  free  surface. 

An  integral  equation  for  the  source  strength  may  be  found  by 
differentiating  (6)  with  respect  to  the  normal  on  the  body  and  setting  it 
equal  to  the  body  boundary  condition  which  requires  that  the  normal  velocity 
on  the  hull  be  zero.  Using  (1)  and  (6),  we  may  write 
1  1  1 

Vnx  =  -  -  a - //  Gn(P,Q)o(Q)dS(Q) - /  vxGn  ( P,  Q)o  ( Q)dn  (7) 

2  4ir  S^j  4rtk0 

Solution  of  (7)  is  achieved  by  discretizing  the  surface  of  the  body  into 
plane  quadrilateral  panels  using  the  method  of  Hess  and  Smith,  as  described  in 
[1]  and  [2].  The  equations  are  then  assembled  to  yield  a  system  of  linear 
algebraic  equations  in  terms  of  the  unknown  source  strength  a  .  The  source 
strengths  and  the  body  configuration  can  then  be  used  as  input  for  the  wave 
amplitude  computations. 

1.2  Determination  of  Wave  Amplitudes 

The  linear  free  surface  elevation  £  at  any  point  (x,y)  is  given  by 
V 

C  =  *  *x  lz  *  0  (8) 

g 

where  the  potential  $  us  given  by  (6).  In  order  to  compute  <fcx  we  need  the 
x-component  of  the  gradient  of  the  Green  function  (4),  which  may  be  decomposed 
as 

(R)  (R')  (P)  (L) 

Gx  *  <JX  +  Gx  +  Gx  +  Gx  .  (9) 

The  first  two  terms  in  (9)  reflect  the  influence  of  the  1 /r  and  1/r'  terms 
(the  Rankine  terms)  in  (4),  and  cancel  out  on  the  z  =  0  plane.  Therefore, 
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we  need  to  evaluate  only  the  terms  due  to  the  double  integral  over  the  body 
surface  (G^PM  and  the  terms  due  to  the  single  integral  along  the 

waterline  Ch(G^lM.  Sy  representing  each  panel  by  a  single  concentrated 
source  at  its  centroid  (xc,yc,zc)  we  may  find 
2 

(P)  2koA  */2 

Gx  - -  /  sec  ?  Imf  FQ  [kQ  sec?Q(ZF  +  Zc)Pd3  ,  (10) 

v  -7i  /  2 

2 

(L)  2kovx^yL  */2 

Gx  - -  r  sec3f'Im{F0(k0  sec?r(ZF  +  ZL)]}dH  .  (11) 

77  -If  /  Z 

Fo(z)  is  given  by 


F0(z)  =  ez  {E-|(z)  -  2ti  i  H[ Im( -z) ] ' 


(12) 


where  H  is  the  Heaviside  step  function  and  E-|(z)  the  complex  exponential 
inteoral 

o  e~ - 

E-,  (z)  =  '  —  dt  .  (13) 

z  t 


In  equations 

(10)  and  (ID, 

A  is 

the  area  of  the  source  panel,  and 

N 

TJ 

II 

i(xF  cosi 

+  yF 

sin-' ) 

9 

(14) 

Zc  =  zc 

-i(xc  cost 

+  yc 

sin'- ) 

9 

(15) 

ZL  = 

-i ( xl  cos? 

+  YL 

sine  ) 

9 

(16) 

(xF,yF,0)  =  the  chosen  field  point  on  z  =  0 

(xc,yc,zc)  =  the  centroid  of  the  source  panel 

<xL»yL*°)  =  the  centroid  of  the  line  segment  of  the  source  panel  that 

coincides  with  the  free  surface 

t.y^  =  the  transverse  component  of  this  line  segment 

vx  =  the  x-component  of  the  unit  normal  of  the  segment. 

The  sign  correction  for  AyL  is  derived  from  the  fact  that  the  line  integral 

in  (7)  is  performed  in  a  counter-clockwise  sense. 


II.  NUMERICAL  CONSIDERATIONS 

The  majority  of  effort  in  the  numerical  evaluation  of  the  wave  amplitudes 
is  spent  in  the  computation  of  the  9-integrals  in  equations  (10)  and  (11). 

To  facilitate  the  computation  we  denote  the  0-integral  by 
it/2 

I  =  -  /  sec30 Im{ F0{ [z+i (x  cos©  +  y  sin0 ) ]sec20} }d0  ,  (17) 

-ir/2 

where  (x,y,z)  is  the  dimensionless  ( non-dimensionalized  by  kQ  )  distance 
between  source  and  field  point.  The  main  difficulty  in  evaluating  (17)  is  due 
to  the  behavior  of  the  Fo<0)  function  versus  0  .  Typical  graphs  of  Fo(0  ) 
are  shown  in  Figures  1  through  6,  and  the  following  comments  can  be  made: 

(i)  F0( 0 )  is  a  rapidly  oscillating  function  in  0  when  z  is  close 

to  zero  (i.e.,  the  source  point  approaches  the  free  surface),  or 
when  x  and  y  are  (in  absolute  value)  large  compared  to  z  . 

This  shows  that  it  is  more  difficult  to  evaluate  the  far-field 
wave  pattern  and  especially  the  waves  generated  by  the  panels 
located  near  to  the  free  surface,  which  give  the  largest 
contribution  to  the  radiated  waves.  This  fact  renders  wave- 
amplitude  computations  more  sensitive  to  the  0 -integration  (17) 
than  the  wave-resistance  (Neuman-Kelvin)  problem. 

(ii)  Fo(0)  is  an  even  function  in  9  ,  i.e.,  Fo(-0)  =  Fo(0)  ,  when 
y  =  0  .  Even  off  the  centerline,  y  =  0  ,  the  graph  of  Fo(0) 
vs.  0  is  almost  symmetric  with  respect  to  0=0  when 

|  x  |  >>  | y |  ,  i.e.,  far  downstream. 

(iii)  Fo(0)  is  an  odd  function  of  y  ,  which  means  that  the 
contribution  from  the  line  integral  (11)  is  zero  right  on  the 
ship's  centerline.  This  fact  is  in  agreement  with  symmetry 
considerations.  Off  the  centerline,  the  line  integral 
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Figure  2.  I»{Fo(0)}  for  (x,y,z)  =  (-10,  0,  0) 


Figure  4.  Im{F0(e)i  for  (x,y,z)  =  (-10,  5,  -0.05) 
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Figure  6.  Im'F0(A)i  for  (x,y,z)  =  (2,  5,  0) 


contribution  approaches  zero  as  the  field  ooint  moves  far 
downstream. 

The  integrand  (17)  possesses  two  integrable  singularities  at  ft  =  *  m/2 
and  another  integrable  one  when  z  =  0  at  ft=  -  tan"’-  (x/y)  due  to  the 
logarithmic-type  singularity  of  the  complex  exponential  integral  E-|(0)  .  In 
the  following  numerical  computations,  the  distance  between  source  and  field 
point  is  changed  from  Cartesian  (x,y,z)  to  cylindrical  (R,i!i,Z)  as  shown  in 
Figure  7.  In  addition,  the  limits  of  integration  in  (17)  are  taken  tc  be 
-n/2  +  e  and  m/2  -  c  ,  where  c  is  a  small  cut-off  number.  Values  of 
integral  (17)  versus  the  number  of  points  in  the  ''-integration,  Np  ,  are 
shown  in  Figures  8  and  9  for  different  values  of  r  and  source-field  point 
distance.  From  these  figures  it  can  be  seen  that:  (a)  integral  (17)  is  well 
behaved  in  the  limit  of  e  -*•  0  and  the  result  is  insensitive  to  the  actual 
value  of  £  ;  and  (b)  when  the  field  point  is  located  off  the  x-axis  the 
convergence  of  (17)  is  much  slower. 

Based  on  the  above  remarks  we  can  compute  integral  (17)  for  a  range  of 
(x,y,z)  values  expressed  in  the  cylindrical  coordinate  system  (R,f,z)  of 
Figure  7.  The  computations  have  been  performed  for  the  following  range  of 
(R,'*i,z)  values: 

R  ranges  from  1  to  a  100  with  increments  of  1. 

p  ranges  from  90°  to  -60°  with  increments  of  10°,  and  from  -60°  to  -90° 
(inside  the  Kelvin  cusp  line)  with  increments  of  0.5°. 

z  takes  the  values  -5,  -2,  -1.5,  and  from  -1.0  to  -0.1  with  increments 
of  -0.1 . 

The  resulting  values  of  integral  (17)  (called  "influence  coefficients”) 
have  been  tabulated  and  can  be  used  to  evaluate  the  influence  of  a  source 
panel  on  any  field  point  within  the  above  range  of  (R,*,z)  values.  Typical 
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variation  of  the  influence  coefficients  versus  z  is  shown  in  Figures  10,  11, 
and  12.  Values  at  z  =  0  have  been  estimated  using  parabolic  extrapolation. 
Rapid  decay  of  the  influence  coefficients  versus  R  outside  of  the  Kelvin 
angle  is  demonstrated  in  Figure  13,  whereas  their  persistence  vs.  R  for  6 
in  the  vicinity  of  the  Kelvin  cusp  (“  -  70.72°)  is  shown  in  Figures  14  and 
15.  The  same  trend  is  evident  from  the  graphs  versus  6  in  Figures  16,  17, 
and  18.  Finally,  the  influence  of  any  source  point  on  the  particular  field 
point  is  estimated  by  using  the  above  tabulated  results  and  trilinear 
interpolation  in  the  (R,A,z)  coordinates.  This  technique  is  only  an 
approximate  one  since  interpolation  errors  are  inevitable.  It  provides, 
however,  an  efficient  and  low  cost  computation  of  a  ship's  Kelvin  wake. 


INTEGRAL  VALUE 


Figure  7.  Coordinate  transformation 


Figure  8.  Convergence  of  integral  (17)  for  (R,i',z)  =  (10,  -90e,  -0.1 


INTEGRAL  VALUE 
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Figure  12.  Influence  coefficients  vs 


z  for  R  =  100 


Figure  13.  Influence  coefficients  vs.  R  for  z  =  -0.1 
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Figure  14.  Influence  coefficients  vs 


R  for  z  =  -0.1  , 


A 


-70 


INFLUENCE  COEFFICIENTS 
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Figure  16.  Influence  coefficients  vs 


•*>  for  z  *  -1  ,  R  =  10 


Figure  18.  Influence  coefficients  vs 


A  for 


-1  , 


100 


III.  RESULTS 


Kelvin  wake  calculations  for  the  QUAPAW  (see  Table  1)  at  three  Froude 
numbers  are  presented.  Source  strengths  were  computed  using  the  Neumann- 
Kelvin  solver  DOCTORS. N-K.  The  original  data  file  XYZFS  was  used  to  create 
quadrilateral  panel  data  for  the  QUAPA’,7.  This  resulted  in  192  panels  per  ship 
side;  24  in  the  longitudinal  direction  and  8  in  the  vertical.  In  the 
following  pages  contour  plots  and  longitudinal  wave  cuts  are  included. 
Computations  performed  by  WAVEAMP  required  less  than  2  seconds  CPU  time  per 
field  point  on  an  Apollo  DN  3000  mini  computer. 

Table  I 

QUAPAW  Model  Dimensions 


LWL  =  16.250  ft 

Wetted  Surface  *  64.883  ft’ 

Beam  *  3.208  ft 

Scale  Ratio  =  12.00 

Draft  -  1 .187  ft 

'VA'VtZ 


Figure  V).  QtlAPAW  near-field  wave  pattern,  Fn  =  0.32 


MODEL  3531  Fn=0.32 
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Wave  cut 


QUAPAW 
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Figure  22.  QUAPAW  near-field  wave  pattern, 
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Wave  cut 


Figure  24.  QUAPAW 


TMT/*Z 


QUAPAW  near-field  wave 
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TMT/*3 


Fiqure  26.  QUAPAW  far-field  wave 


IV.  PROGRAM  DOCUMENTATION 


IV. 1  Description  of  Input 

WAVEAMP  needs  three  data  files  to  run.  File  #1  contains  the  influence 
coefficients  and  is  read  through  unit  1.  File  #2  contains  wave  resistance 
data,  hull  and  panel  data,  and  source  strengths.  It  is  read  in  on  unit  5. 

This  file  is  produced  by  running  the  Heumann-Kelvin  solver  DOCTORS. N-K.  The 
user  needs  to  prepare  one  data  file.  File  #3,  to  be  read  on  unit  4.  This  file 
contains  the  coordinates  of  the  field  points  and  is  read  through  free  format 
statements.  The  first  line  contains  the  total  number  of  field  points.  There 
is  no  limitation  on  the  maximum  number  of  points.  Each  subsequent  line 
contains  the  (x,y,z)  coordinates  of  the  particular  field  point  (in  feet)  in 
that  order.  The  following  sign  convention  applies:  x  is  zero  at  the  ship's 
mid-section,  positive  upstream,  and  negative  downstream;  y  is  positive 
starboard;  and  z  is  zero. 

IV. 2  Description  of  Output 

In  addition  to  on-screen  information  output,  WAVEAMP  produces  one  results 
file  written  on  unit  6.  This  file  is  identical  to  input  File  #3  except  that 
the  third  column  z  is  substituted  by  c(x,y)  ,  the  wave  elevation  (in  feet) 
at  the  particular  field  point.  The  wave  elevation  is  positive  for  a  wave 
crest  and  negative  for  a  wave  trough. 

IV. 3  Run  Commands 

The  program  is  written  in  standard  Fortran-77  for  an  Apollo  mini-computer 
and  can  be  easily  transferred  to  any  Fortran  supporting  machine.  Execution 
begins  (for  the  Apollo  machine)  by  simply  typing  the  name  of  the  executable 
code  and  it  prompts  the  user  for  the  data  and  results  files  names. 
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IV. 4  Error  Messages 

The  only  case  where  error  messages  may  appear  when  running  WAVEAMP  is 
when  the  distance  between  source  and  field  point  is  such  that  (R,4,z)  are 
out  of  their  tabulated  values  range,  as  described  in  Section  II.  In  such  a 
case,  execution  is  terminated  and  a  self-explanatory  error  message  is  printed 
at  the  end  of  the  results  file. 


CONCLUSIONS 


The  procedure  of  tabulating  the  values  of  the  influence  coefficients 
adopted  by  program  WAVEAMP  is  an  efficient  and  inexpensive  alternative  to  time 
consuming  but  more  accurate  computations  of  the  waves  generated  by  a  body 
moving  in  an  inviscid,  calm  fluid.  The  achieved  computational  efficiency 
outweighs  errors  introduced  by  the  interpolation  scheme  and  simple  monopole 
method.  Finally,  computational  accuracy  may  be  improved  by  using  a  finer  grid 
of  tabulated  values  and  more  accurate  interpolation  schemes. 
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APPENDIX:  PROGRAM  LISTING 


The  following  pages  provide  a  complete  listing  of  program  WAVEAMP.  The 
comments  incorporated  in  the  source  code  should  be  adequate  to  describe  the 
flow  of  the  code. 


non  noon  ooo  o  on  o  ooocjooooooooooo 


waveamp . ftn 


06/07/88 


subroutine:  Sinain 


PROGRAM  WAVEAMP 

VERSION  3  :  JUNE  1988 

PROGRAMMER  :  F . A.  PAPOULIAS 

DATE  :  FEBRUARY  1988 

LOCATION  :  THE  UNIVERSITY  OF  MICHIGAN 

COMPILER  :  STANDARD  FORTRAN-77 

APOLLO  COMPUTER 


CALCULATION  OF  WAVE  ELEVATIONS  CREATED  BY  A  THICK  SHIP 
TRAVELLING  AT  A  CONSTANT  SPEED  BY  THE  NEUMANN -KELVIN  METHOD 
MONOPOLES  APPROXIMATION 

DIMENSION  VEE (200) ,SIGMA(200) ,XXO(3, 4,200) , 

1  XXQ (3,4,200)  ,  XXC  ( 3 ,  200)  ,  TRANS < 9,  200 ) , ANN (3, 200) , 

2  D  (4, 200) , TMAX (200)  ,XXE(3,  4,200),AA(8,200),AM(4,200) 

3  XXF (3,200)  ,  ICNT (4)  ,  WVCT (3) , XXP (3,4,200) 

COMMON/DATA1/X1A (14)  ,  X2A ( 100 ) , X3A (7  6) 

COMMON/DATA2/YA (14,  100,  7  6) 

DATA  PI/  3.1415  92653  5898/ 

OPEN  (UNIT-1,  FILE-' ‘INFLUENCE  COEFF .  FILE-') 

OPEN  (UNIT-4,  FILE-' ‘FIELD  POINTS  FILE-') 

OPEN  (UNIT-5,  FILE-' ‘WAVE  RES.  DATA  FILE-') 

OPEN  (UNIT-6,  FILE-' ‘RESULTS  FILE-') 

WRITE  (*,1002) 

READ  WAVE  RESISTANCE  DATA  FILE 


READ  (5,3001)  I 
READ  (5,3004) 
READ  (5,3002) 
READ  (5,3002) 
READ  (5,3002) 
READ  (5,3002) 
READ  (5,3002) 
READ  (5,3002) 
READ  (5,3002) 
READ  (5,3002) 
READ  (5,3002) 
READ  (5,3002) 
READ  (5,3002) 
READ  (5,3003)  i 
READ  (5,3002) 
READ  (5,3002) 

WRITE  (*,1003) 


NP , NX, NZ, NV,  NH,  IFREE, IPANEL, ISP,  ISL 
TOLL 

(SIGMA(IP) , IP-1, NP) 

(  (  (XX0(I3, 14,  IP)  ,  IP-1,  NP)  ,  14=1, 4)  ,  13=1, 3) 
(  (  (XXQ  (13,  14,  IP)  ,  IP-1,  NP)  ,  14  =  1,  4)  ,  13=1, 3) 
( (XXC  (13,  IP)  ,  IP-1,  NP)  ,  13-1,  3) 

( (TRANS  (19,  IP)  ,  IP-1,  NP)  ,  19=1,  9) 

(  (ANN (13,  IP)  ,  IP-1 , NP ) , 13-1,3) 

(  (D(I4, IP)  ,  IP-1,  NP) ,14-1,4) 

(TMAX(IP)  ,  IP-1, NP) 

( (AA(I8, IP)  ,  IP-1, NP) ,18-1,8) 

( (AM (14,  IP),  IP-1,  NP) ,14-1,4) 

(  (XXF (13,  IP)  ,  IP-1 ,  NP ) ,13=1,3) 

RATIOl, RATI02, G, V 

(  (  ( XXE  (13,  14,  IP)  ,  IP-1,  NP)  ,  14  =  1, 4)  ,  13=1, 3) 
(  (  (XXP  (13,  14,  IP)  ,  IP-1,  NP)  ,  14  =  1,  4)  ,  13  =  1,  3) 


READ  INFLUENCE  COEFFICIENTS  DATA  FILE 

DO  1  11-1,14 
DO  2  12-1, 100 

READ  (1,301)  (YA(I1,  12, 176) ,  176  =  1,  76) 
2  CONTINUE 
1  CONTINUE 

LOAD  VECTOR  OF  DEPTH  COORDINATE  Z 

DO  3  1-1,14 

IF  (I.EQ.l)  2— 5.0 
IF  (I.EQ.2)  Z— 2.0 


nno  o  ooo  non  n  non  n  non  non  non  non  n  non  nnn 
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IF  (I.EQ.3)  Z— 1.5 

IF  <  I .  GT .  3 )  Z— 1.0+U-4)  *0.1 

X1A(I)-Z 

3  CONTINUE 

LOAD  VECTOR  OF  RADIAL  COORDINATE  R 

DO  4  1-1,100 
R-I 

X2A ( I ) -I 

4  CONTINUE 

LOAD  VECTOR  OF  ANGULAR  COORDINATE  PHI 

DO  5  1-1,76 

PHI-90. 0-10. 0MI-1) 

IF  (I.GE.16)  PHI— 60. 0-0 . 5*  (1-16) 

X3A ( I ) —PHI 

5  CONTINUE 

READ  (4,  ‘)  NFP 

LOOP  OVER  FIELD  POINTS 

DO  800  KP— 1, NFP 

WRITE  (*,1001)  KP , NFP 
READ  (4,*)  (WVCT(I) , 1-1, 3) 

INITIALISE  VEE  VECTOR 

DO  300  1-1 , NP 
VEE (I) -0 . 0 

300  CONTINUE 

CONSIDER  BOTH  SIDES  OF  THE  SYMMETRY  PLANE  OF  THE  SHIP 
DO  450  IY-1,2 

IF  (IY.EQ.2)  WVCT (2 ) — -WVCT  (2) 

THE  LINE-INTEGRAL  CONTRIBUTION  OF  THE  WAVE  TERM 
FACLIN-V**2/ (4 . 0*PI*G) 

CALL  LINELV  (NP,  NX,  4,  WVCT,  XXQ,  ANN,  G,  V,  FACLIN,  VEE,  IY) 
THE  SURFACE- INTEGRAL  CONTRIBUTION  OF  THE  WAVE  TERM 
FACPAN— 1.0/  (4 . 0*PI) 

CALL  PANELV  (NP,  XXF,  WVCT,  AA,  G,  V,  FACPAN,  VEE) 

IF  (IY.EQ.2)  WVCT (2) --WVCT  (2) 

450  CONTINUE 

COMPUTE  WAVE  ELEVATION 

DO  780  JP— 1 , NP 

WVCT  (3)- (VEE  (JP)  ‘SIGMA  ( JP)  )  *  (V/G) +WVCT  (3) 

"80  CONTINUE 

PRINT  OUT  RESULTS 

WRITE  (6,2000)  WVCT (1) , WVCT (2) , WVCT (3) 

=00  CONTINUE 

STOP 

FORMAT  STATEMENTS 

301  FORMAT  (5E15.5) 

1301  FORMAT  ('  FIELD  POINT', 16,'  OUT  OF' , 16, '  TOTAL') 
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1002  FORMAT  ('  READING  WAVE  RESISTANCE  DATA  FILE 

1003  FORMAT  ('  READING  INFLUENCE  COEFFICIENTS 
2000  FORMAT  (3E15.5) 

3001  FORMAT  (915) 

3002  FORMAT  (5E20.10) 

3003  FORMAT  (4E20.10) 

3004  FORMAT  (E20.10) 

C 

END 

C 
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waveamp.ftn  06/07/88  subroutine:  panelv 


SUBROUTINE  PANELV  (NP, XXC,  WVCT, AAA, G, V, FACPAN, VEE) 

C 

C  THE  WAVE  PART  OF  THE  VELOCITY  INDUCED  AT  A  FIELD  POINT 

C  BY  A  UNIFORM  SOURCE  PANEL 

C 

DIMENSION  XXC (3, NP)  ,  WVCT (3)  ,  AAA (8, NP) , VEE (NP) 

C 

COMMON/DATA1/X1A (14) ,  X2A (100 ) , X3A (7  6) 

COMMON/DATA2/YA(14, 100,76) 

C 

DATA  PI/  3.1415  92653  5898/ 

C 

AK0»G/V**2 

C 

C  STEP  THROUGH  THE  SOURCE  PANELS 

C 

DO  300  JP“1,NP 

DX-WVCT ( 1 ) -XXC ( 1 ,  JP ) 

DY«WVCT<2) -XXC (2, JP) 

DZ-WVCT ( 3 ) +XXC ( 3 , JP ) 

DY-ABS (DY) 

DZ— ABS  (DZ) 

DX-DX*AK0 

DY-DY*AK0 

DZ-DZ*AK0 

IF  (DZ.LT. (-5.0) )  GO  TO  8001 
IF  (DZ .GE .0.0)  GO  TO  8002 

R-SQRT (DX*DX+DY*DY) 

IF  (R . LT .1.0)  GO  TO  8003 
IF  (R.GE. 100.0)  GO  TO  8004 
PHI«ASIN(ABS (DX) /R) 

PHI«PHI*180 .0/PI 
IF  (DX .  LT  .0.0)  PHI— PHI 
CALL  LOCATE (XI A,  14,  DZ,I) 

CALL  LOCATE (X2 A, 100,  R,J) 

CALL  LOCATE (X3A,  76, PHI, K) 

CALL  TRLINT (14,  100 , 7 6,  I ,  J,  K,  YA,  XlA,  X2A, X3A, DZ, R, PHI, SUMX) 
FAC— 2.0*AK0**2*AAA<1,  JP)  /PI 
VXX-FAC* SUMX 

VEE (JP)-VEE (JP) +FACPAN*VXX 
300  CONTINUE 
C 

RETURN 

C 

C  AB-NORMAL  OUTPUT 

C 

8001  WRITE  (6,7001)  JP 
STOP  9001 

8002  WRITE  (6,7002)  JP 
STOP  9002 

8003  WRITE  (6,7003)  JP 
STOP  9003 

8004  WRITE  (6,7004)  JP 
STOP  9004 

7001  FORMAT  (/,'  Z  IS  LESS  THAN  -5  FOR  PANEL  NO. ',15) 

7002  FORMAT  (/,'  Z  IS  GREATER  THAN  0  FOR  PANEL  NO.',  15) 

7003  FORMAT  (/,'  R  IS  LESS  THAN  1  FOR  PANEL  NO. ',15) 

7004  FORMAT  (/,'  R  IS  GREATER  THAN  100  FOR  PANEL  NO. ',15) 

C 

END 

C 
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waveamp. ftn 


06/07/88 


subroutine:  linelv 


SUBROUTINE  LINELV  (NP,  NX,  NV,  WVCT,  XXQ,  ANN,  G,  V,  FACLIN,  VEE,  IFAP) 

THE  WAVE  PART  OF  THE  VELOCITY  INDUCED  AT  A  FIELD  POINT 
BY  A  UNIFORM  SOURCE  LINE 

DIMENSION  XXQ(3,NV,NP)  ,WVCT(3) ,ANN(3,NP) ,VEE(NP) 

C0MM0N/DATA1/X1A(14)  ,X2A(100) ,X3A(76) 

COMMON/DATA2 / YA (14,100,76) 

DATA  PI/  3.1415  92653  5898/ 


AK0=G/V**2 
TOL-1.0E-05 
NX1-NX-1 

STEP  THROUGH  THE  SOURCE  LINES 
DO  300  JX-1,NX1 

DX-WVCT(l) -0.5* (XXQ (1,  4,  JX) +XXQ(1, 1, JX) ) 

DY«WVCT(2) -0.5* (XXQ  (2,  4, JX)+XXQ(2, 1, JX) > 

DZ-WVCT(3) 

DY-ABS (DY) 

DZ— ABS  (DZ) 

DX-DX*AK0 
DY«DY*AK0 
DZ«DZ*AK0 

DDY«XXQ(2, 1, JX) -XXQ (2, 4, JX) 

IF  (ABS (DDY) .LT.TOL)  GO  TO  300 
IF  (DZ . LT . (-5.0))  GO  TO  8001 
IF  (DZ .GE .0.0)  DZ-0.001 
R-SQRT (DX*DX+DY*DY) 

IF  (R.LT.l.)  GO  TO  8003 

IF  (R.GE. 100.0)  GO  TO  8004 
PHI-ASIN(ABS(DX)/R) 

PHI-PHI*180.0/PI 
IF  (DX.LT.0.0)  PHI— PHI 
CALL  LOCATE (X1A,  14,  DZ,I) 

CALL  LOCATE (X2A, 100,  R,J) 

CALL  LOCATE (X3A,  7  6,  PHI, K) 

CALL  TRLINT (14,  100,7  6,  I,  J,  K,  YA, X1A, X2A, X3A, DZ, R, PHI, SUMX) 
FAC-2.0*DDY*AK0**2*ANN<1, JX) /PI 
VX-FAC*SUMX 
IF  (IFAP.EQ.2)  VX— VX 
VEE ( JX) -VEE ( JX) +FACLIN*VX 
0  CONTINUE 

RETURN 

AB-NORMAL  OUTPUT 

e001  WRITE  (6,7001)  JP 
STOP  9001 

8003  WRITE  (6,7003)  JP 
STOP  9003 

8004  WRITE  (6,7004)  JP 
STOP  9004 

7001  FORMAT  (/,'  Z  IS  LESS  THAN  -5  FOR  PANEL  NO. ',15) 

7003  FORMAT  (/,'  R  IS  LESS  THAN  1  FOR  PANEL  NO. ',15) 

7004  FORMAT  (/,'  R  IS  GREATER  THAN  100  FOR  PANEL  NO. ',15) 

C 

END 

C 
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waveamp.ftn  06/07/88  subroutine:  locate 


SUBROUTINE  LOCATE  (XX,N,X,J) 

GIVEN  AN  ARRAY  XX  OF  LENGTH  N,  AND  GIVEN  A  VALUE  X,  RETURNS 
A  VALUE  J  SUCH  THAT  X  IS  BETWEEN  XX (J)  AND  XX(J+1)  .  XX  MUST 
BE  MONOTONIC,  EITHER  DECREASING  OR  INCREASING.  J=0  OR  J=N 
IS  RETURNED  TO  INDICATE  THAT  X  IS  OUT  OF  RANGE 

DIMENSION  XX (N) 

C 

JL-0 

JU-N+1 

10  IF  (JU-JL.GT . 1)  THEN 
JM= (JU+JL) /2 

IF  ( (XX(N) .GT.XX(l) ) .EQV. (X.GT.XX(JM) ) )  THEN 
JL=JM 
ELSE 
JU=JM 
ENDIF 
GO  TO  10 
ENDIF 
J=JL 
RETURN 
END 
C 


ooooooooo 
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waveamp.ftn  06/07/88  subroutine:  trlint 


SUBROUTINE  TRLINT (L,M,  N,  I,  J,  K, YA, X1A,  X2A, X3A, XI,  X2, X3,  Y> 

LINEAR  INTERPOLATION  IN  THREE  DIMENSIONS. 

GIVEN  A  THIRD  ORDER  TENSOR  YA  (L, M, N) ;  ARRAYS  X1A(L), 

X2A (M) ,  X3A (N)  ;  THREE  NUMBERS  X1,X2,X3  AND  I, J, K  SUCH 
THAT  XI  IS  BETWEEN  X1A<I)  AND  X1A(I  +  1),  X2  IS  BETWEEN 
X2A  ( J)  AND  X2A  ( J+l)  ,  AND  X3  IS  BETWEEN  X3A(K)  AND 
X3A (K+l)  ;  IT  RETURNS  Y(X1,X2,X3)  BY  TRILINEAR 
INTERPOLATION  IN  THE  TENSOR  YA. 

DIMENSION  YA{L,M,N)  ,  X1A  (L)  ,  X2A  (M)  ,  X3A  (N) 

C 

Y1-YA(I  ,J  ,  K  ) 

Y2-YA(I  , J+l , K  ) 

Y3=YA(I  , J+l, K+l) 

Y4«YA(I  ,J  ,K+1) 

Y5-YA(I+1,J  , K  ) 

Y6=YA ( 1+1, J+l,  K  ) 

Y7=YA (1+1, J+l, K+l) 

Y8«=YA(I+1,J  , K+l ) 

C 

T* (XI -XI A ( I ) ) / (XI A (1+1)  -X1A(I)  ) 

U“ (X2-X2A ( J) ) / (X2A ( J+l ) -X2A ( J)  ) 
W*(X3-X3A(K))/(X3A(K+1)-X3A(K)  ) 

C 

Tl=l . 0-T 
U1-1.0-U 
W1-1.0-W 

C 

Y=T1*U1*W1*Y1+T*U1*W1*Y2+T*U*W1*Y3+T1*U*W1*Y4  + 

.  T1*U1*W  *Y5+T*U1*W  *Y6+T*U*W  *Y7+T1*U*W  *Y8 

C 

RETURN 

END 


